#

push!(LOAD_PATH, "/home/mary/workspace/")
using pimary

#using Plots
include("../test/visliz_helper.jl")

using LinearAlgebra
BLAS.set_num_threads(1)

#=
rs = a / aB
#是平均距离
#在球形的区域 a^3 4π/3 = V/N -> a = (3/4πn)^(1/3)
#在正方形的区域 a^3 8 = V/N -> a = 1/2 (1/n)^(1/3)

假设rs = 2 ，将aB取成1则a=2，想包含64个电子进去的话，那么总的面积就是
64 = V/N
V = 64*64
正方形区域的直径就是
R = V^{1/3} = 16

rs = 1/2 (1/n)^{1/3}
n =  1 / (2rs)^3

#温度 Θ = kB T/ Ef
Kf^3 / (1/R)^3 = Np
Kf = Np^{1/3} / R
Ef = Kf^2 / 2m
#
1/β = kB T = Θ*Ef = Θ * (Np^{1/3} / R)^2 / 2m
=#

#1bohr=0.529177Å = 0.529177e-8cm
#M是摩尔质量（g/mol），Na阿佛加德罗常数（6.022*1023/mol），ρ是质量密度（g/cm3）
#M=1, ρ=1.225 kg/m3 = 1225 g/ 1000000 cm3 = 0.001225 g/cm3
#则原子数密度为，na=Na*ρ/M(cm-3)=0.08923673ρ/M(bohr-3)
#na = 0.00010931499425
#L=(N/na)1/3

"""正方"""
function rs_to_R(Np)
    #rs = 1.86 对应 na = 2.503e23 cm-3 = 2.503e23 (10^-8)^3 Å-3
    #na = 0.2503 Å-3 = 0.2503 * (0.529177)^3 bohr
    #na = 0.03709058893420576 bohr-3
    na = 0.03709058893420576
    R = (Np/na)^(1/3)
    return R
end


"""从Θ到beta"""
function k_to_β(k)
    #1K相应于3.16679*10-6Eh
    Eh = k*3.16679e-6
    return 1/Eh
end



function worldlines(betaval)
    rs = 2
    Np = 9
    Nc = 2*Np
    Θ = 8
    Nt = 10
    effR = rs_to_R(Nc)
    totVol = effR^3
    Dprep = (totVol/Nc)^(1/3)
    println("effR = $(effR); D = $(Dprep)")
    Xpos = [3, 3, 2]
    @assert prod(Xpos) == Nc
    R = Xpos * Dprep
    println("R = $(R)")
    β = k_to_β(betaval)
    Δβ = β/Nt
    println("β = $(β), Δβ = $(Δβ)")
    #
    λe = sqrt(2π*(pimary._hbar^2)*Δβ/pimary._ElecMass)
    λ11 = sqrt((pimary._hbar^2)*Δβ/pimary._ElecMass)
    λ22 = sqrt((pimary._hbar^2)*Δβ*(1/pimary._ElecMass+1/pimary._NeurMass)/2)
    println("lambda_e=$(λe) lambda_11=$(λ11) lambda_22=$(λ22)")
    #绘制当前势能的图像
    xvals = collect(Base.OneTo(100))/25
    pote2 = [Kelbg_attract([x, 0.0, 0.0], [2.0, 0.0, 0.0]; λ22=λ22, Cne=1.0) for x in xvals]
    #pltp = plot(xvals, pote2, xlims=(0.0, 4.0), ylims=(-20.0, 0.0))
    #savefig(pltp, "potential_jelluim.png")
    #生成初始的位置
    R_c = zeros(Nc, 3)
    cidx = 0
    for x in 1:1:Xpos[1]; for y in 1:1:Xpos[2]; for z in 1:1:Xpos[3]
        cidx += 1
        R_c[cidx, 1]  = (x-1) * Dprep
        R_c[cidx, 2]  = (y-1) * Dprep
        R_c[cidx, 3]  = (z-1) * Dprep
    end; end; end
    fbne = WrdFb(Float64, Nc, Nt)
    for ip in Base.OneTo(Nc); for it in Base.OneTo(Nt)
        fbne[ip][it, :] = map_to_1st(R_c[ip, :], R)
    end; end
    #plt = plot()
    #draw_fb3D!(plt, fbne, :black)
    #随机的电子位置
    fbup = random_path(Float64, Np, Nt, R; λe=λe)
    fbdn = random_path(Float64, Np, Nt, R; λe=λe)
    #draw_fb3D!(plt, fbup, :red)
    #draw_fb3D!(plt, fbdn, :blue)
    #savefig(plt, "test_jellium.png")
    #
    rhov, rhop, rhok = spinful_rho_path(
        fbup, fbdn, fbne, R; λe=λe, μee=λ11, μep=λ22, Cne=1.0
    )
    println(rhov)
    println(rhop)
    println(rhok)
    #
    for it in Base.OneTo(Nt)
        upbds1 = beads(fbup, it)
        nebds1 = beads(fbne, it)
        intmat = pimary.Kelbg_attract_pairs2(upbds1, nebds1; λ22=λ22, Cne=1.0)
        dis = sqrt(sum((upbds1[1, :]-nebds1[1, :]).^2))
        #println(dis)
        scatter!(pltp, [2-dis], [intmat[1, 1]])
    end
    savefig(pltp, "potential_jelluim.png")
    #开始随机的抽样过程
    #
    nbin = 10
    nmeas::Int = 1e5
    #WARM UP
    for _ in Base.OneTo(nmeas)
        pimary.update!(fbup, fbdn, fbne, R; λe=λe, μee=λ11, μep=λ22, Cne=1.0)
    end
    #
    sgn_tot = zeros(nbin)
    eng_tot = zeros(nbin)
    for bidx in Base.OneTo(nbin)
        #plt = plot()
        #draw_fb3D!(plt, fbne, :black)
        #draw_fb3D!(plt, fbup, :red)
        #draw_fb3D!(plt, fbdn, :blue)
        #savefig(plt, "test_$(bidx)_jellium0.png")
        #
        for idx in Base.OneTo(Int(nmeas))
            pimary.update!(fbup, fbdn, fbne, R; λe=λe, μee=λ11, μep=λ22, Cne=1.0)
            engv = pimary.internal_energy(fbup, fbdn, fbne, R, Δβ; λe=λe, μee=λ11, μep=λ22, Cne=1.0)
            eng = engv[1] + engv[2] + engv[3] - engv[4]
            eng += 1.5 * (Nc + Np*2) / β
            #println(engv, eng)
            sgn = pimary.abs_phase(fbup, fbdn, R; λe=λe)
            #println(sgn)
            sgn_tot[bidx] += sgn
            eng_tot[bidx] += sgn*eng
            #_, _, rk = spinful_rho_path(fbup, fbdn, fbne, R; λe=λe, μee=λ11, μep=λ22, Cne=1.0)
            #println("rk $(idx): ", sum(rk))
            #if mod(idx, 1e3) == 0
            #    plt = plot()
            #    draw_fb3D!(plt, fbne, :black)
            #    draw_fb3D!(plt, fbup, :red)
            #    draw_fb3D!(plt, fbdn, :blue)
            #    savefig(plt, "test_$(bidx)_jellium$(idx).png")
            #    println(sgn_tot)
            #    println(eng_tot)
            #end
        end
    end
    #
    for bidx in Base.OneTo(nbin)
        eng_tot[bidx] = eng_tot[bidx] / sgn_tot[bidx]
        sgn_tot[bidx] = sgn_tot[bidx] / nmeas
    end
    println(sgn_tot)
    println(eng_tot)
    aveeng = sum(eng_tot)/nbin
    erreng = @. (eng_tot - aveeng)^2
    erreng = sqrt(sum(erreng)/(nbin-1))
    println(aveeng)
    println(erreng)
end



#println(rs_to_R(2.0, 64))
#println(Θ_to_β(8, 2.0, 64, 1.0))
worldlines(1e4)

